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Applications  and  time-domain  solution  of  higher-order  parabolic 
equations  in  underwater  acoustics 

Michael  D.  Collins31 

Naval  Ocean  Research  and  Development  Activity,  Stennis  Space  Center,  Mississippi  39529 
( Received  9  February  1989;  accepted  for  publication  19  May  1989) 

A  higher-order  parabolic  equation  and  the  corresponding  higher-order  time-domain  parabolic 
equation  are  derived  from  a  Pade  series  and  solved  numerically.  These  models  provide  accurate 
solutions  for  problems  involving  very-wide-angle  propagation  (e  g.,  propagation  in  the 
nearfield  or  over  a  hard  ocean  bottom);  propagation  in  domains  in  which  sound-speed 
variations  are  large  (e.g.,  propagation  in  deep  water,  deep  within  the  ocean  bottom,  in  high¬ 
speed  ocean  bottoms,  and  possibly  ,,f different  wave  types);  and  propagation  out  to  very  long 
ranges.  The  possibility  of  modeling  elastic  wave  propagation  with  a  similar  approach  is 
considered. 

PACS  numbers:  43.30.Bp,  43.20.Bi 


INTRODUCTION 

The  parabolic  eo.uation  (TE)  has  undergone  extensive 
development  since  it  was  first  applied  to  underwater  acous¬ 
tics'  including  improvements  in  accuracy  and  implementa¬ 
tion  in  the  time  domain.  The  phase  errors  of  PE  solutions, 
which  approximate  the  solution  of  the  wave  equation,  were 
reduced  significantly  with  the  introduction  of  the  wide-angle 
PE.2-4  Although  various  generalizations  of  the  wide-angle 
PE  have  been  considered, 5-7  the  aperture  limitation  of  the 
PE  has  remained  an  issue  of  concern.  The  time-domain  para¬ 
bolic  equation  (TDPE)8'12  allows  one  to  perform  pulse 
propagation  calculations  without  Fourier  synthesis.  The 
TDPE  has  been  extended  to  handle  interface  conditions,9 
nonlinear  propagation, 10  density  variations  and  sediment  at¬ 
tenuation,1 1  and  wide-angle  diffraction  and  sediment  disper¬ 
sion.12 

In  this  article,  a  higher-order  PE  based  on  a  Pade  series7 
is  shown  to  provide  solutions  comparable  in  accuracy  to  nor¬ 
mal-mode  solutions  for  problems  involving  very-wide-angle 
propagation,  large  variations  in  sound  speed,  and  propaga¬ 
tion  out  to  long  range.  Since  the  Pade  series  is  composed  of 
rational  linear  terms,  the  higher-order  PE  is  easy  to  solve 
numerically.  The  corresponding  higher-order  TDPE  is 
solved  numerically  and  compared  with  a  wide-angle  TDPE 
designed  for  propagation  in  shallow  water. 12  The  possibility 
of  applying  the  Pade  series  to  derive  a  PE  for  elastic  wave 
propagation  is  investigated. 

I.  THE  HIGHER-ORDER  PE 

A  time-harmonic  steady  state  is  assumed,  and  the 
acoustic  pressure/?  is  factored  asp(x,/)  =  P(x)ex p(  —  icot), 
where  t  is  time,  x  is  the  Cartesian  position  vector,  and  to  is  the 
circular  frequency.  The  complex  pressure  P  is  assumed  to 
satisfy  the  pressure-release  boundary  condition  P  =  0  at  the 
ocean  surface,  the  outgoing  radiation  condition  at  infinity, 
and  the  reduced  wave  equation 

pV-[(\/p)VP]  +  K:P=  —  4n"<5(x  —  x„),  (1) 
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where  the  point  x()  is  the  source  location.  The  complex  wave- 
number  K  —  k  +  i<tP  j  is  used  to  account  for  sediment  lose. 
The  wavenumber  is  k  =  co/c,  t]  =  (40rr  log,,,  e)  ',p  is  the 
attenuation  in  decibels  per  wavelength  (dB/A  ),p  is  the  den¬ 
sity,  and  c  is  the  sound  speed. 

We  assume  that  variations  in  azimuth  are  negligible  and 
solve  Eq.  ( 1 )  in  cylindrical  coordinates,  with  z  being  the 
depth  below  the  ocean  surface  and  r  being  the  horizontal 
distance  from  a  source  at  the  depth  z„.  Variations  in  range 
are  assumed  to  be  sufficiently  weak  so  that  dp/dr  can  be 
ignored,  which  simplifies  Eq.  ( 1 )  to 

dr  p  3z  dz  dr  r  dr 

=  -  (2/r)8(r)8(z-z0).  (2) 


We  define  Q  =  JrP,  and  for  r>  0  Eq.  (2)  becomes 

dzr  p  dz  dz  dr  4  r 
We  assume  that  r>  r0,  where  kr0 >  1,  and  drop  the  0(r~2) 
term  in  Eq.  (3)  to  obtain  the  farfield  equation 


dz 2  p  dz  dz  dr 


(4) 


For  range-independent  domains,  Eq.  (4)  factors  exactly  to 


where  c„  is  a  reference  sound  speed,  k„  —  co/c„,  and 


dp_d_ 
dz 2  p  dz  dz 


(6) 


Equation  (5),  which  we  refer  to  as  PE,  ,  is  an  accurate  ap¬ 
proximation  for  many  range-dependent  problems  in  under¬ 
water  acoustics  and  can  be  solved  in  terms  of  outgoing  cou¬ 
pled  modes." 

PE  derivations  are  based  on  approximations  of  the  func¬ 
tion  v  1  a-  .r  1  <xr  i  vj  ^  i  <-ich  as  the  Taylor  series 

v'  1  +  Jt  -  1  =  U  -  J*2  +  ,'x'  -  ■  ■  •  .  (7) 
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If  K  =  kt),  the  first  term  of  this  series  and  the  approximation 

(K2  -  *  j;  )/2k„ sk-k„  +  ii/p \k  |  ( 8 ) 

give  the  narrow-angle  PE 

*E  =  Hk-kfl)U-n0\k\U  +  ^U,  (9) 

or  2k, , 

where  Q—U exp(/A',/).  The  terms  on  the  right  side  of  Eq. 
(9)  (from  left  to  right)  are:  the  refraction  term,  which  ac¬ 
counts  for  variations  in  sound  speed;  the  loss  term,  which 
accounts  for  sediment  attenuation;  and  the  diffraction  term, 
which  accounts  for  the  vertical  component  of  propagation. 1,1 

The  first  two  terms  of  the  Taylor  series  have  been  used  to 
derive  a  wide-angle  PE.'  However,  the  higher-order  PE’s 
based  on  the  Taylor  series  are  relatively  inefficient.  Since  the 
Taylor  series  diverges  for  |xj>  1,  many  terms  in  the  series 
are  required  for  problems  involving  large  differences  in 
sound  speed  or  very-wide-angle  propagation.  Furthermore, 
higher-order  PE’s  based  on  the  Taylor  series  are  difficult  to 
implement  because  x  is  raised  to  powers,  which  results  in 
operators  raised  to  powers. 

The  wide-angle  PE  is  based  on  a  idtional  linear  Pade 
approximation.  A  generalization  of  the  wide-angle  PE  that  is 
based  on  a  ratio  of  polynomials'’  is  difficult  to  implement 
because  powers  of  x  are  involved.  A  generalization  of  the 
wide-angle  PE  that  involves  only  first  powers  is  based  on  the 
following  Pade  series7: 

v  1  -r  .v  ~  1  =  V  a"'X  +  0(x2"  ’  ’),  (10) 

i  i  1  +  bIMx 

where  n  is  the  number  of  terms  in  the  Pade  expansion  and 

ajn  —  [2/(2h  +  1 )  ]sin’[  jir/(2n  +  1 )  J,  (11) 

cos:[y-/(2«  +  1)].  (12) 

Since  the  Pade  series  is  valid  outside  the  radius  of  con¬ 
vergence  of  the  Taylor  series,  relatively  few  terms  are  needed 
for  ix|  »  1.  We  illustrate  this  in  Table  I.  The  four-term  Tay¬ 
lor  series  is  better  than  the  one-term  Pade  series  forx  <  1,  but 
the  one-term  Pade  series  is  better  for  x>  1.  The  two-term 
Pade  series  and  the  four-term  Taylor  series  are  both  correct 
to£7(.t')  for  small  x.  Yet  the  two-term  Pade  series  is  substan¬ 
tially  better  than  the  four-term  Taylor  series.  The  three-term 


TABL.F.  I  Comparison  of  Taylor  and  Pade  series. 


X 

Four-term 

Taylor 

One-term 

Fade 

Two-term 

Fade 

1  hree-term 
Pade 

*  1  4  .x  1 

0.2? 

0  11X01 

O  1 17  65 

0  118  03 

0.118  03 

0.118  03 

0.50 

0 224  12 

0.222  22 

0  224  72 

0.224  74 

0.224  74 

0.7S 

O  MS  70 

0  315  74 

0.322  74 

0  322  87 

0.322  88 

1  (XI 

0..VW  44 

0  4001X) 

0  413  79 

0  414  20 

0.414  21 

1.25 

0  450  39 

0  476  19 

0.4 W  04 

0.499  96 

0.500  (X) 

t  .50 

0.4X1 

0.545  45 

0.579  31 

0.581  05 

0  581  14 

1.75 

0.460  7X 

0.608  7i) 

0.655  23 

0.658  12 

0  658  31 

MX) 

0.375 (X) 

0.720  3X) 

0.795  84 

0.802  21 

0.802  78 

2  50 

0.0XO  <7 

o  7  1 

. 

::  ov  v4 

0  870  83 

2. ;  5 

-  w.Xh  ./j 

O.XI4  81 

0.923  76 

0.935  19 

0.936  49 

3.00 

1.101  56 

0.857  14 

0.983  61 

0.998  17 

[.(XXI  rxi 

Pade  series  is  fairly  accurate  well  beyond  the  radius  of  con¬ 
vergence  of  the  Taylor  series  near  x  =  3. 

The  Pade  series  gives  the  higher-order  PE 

=  *„  y  v,  ,m 

dr  i  i  k  +  bj  „  (L  +  K~  —  k;,) 

which  we  refer  to  as  PE„ .  Equation  (13)  can  be  solved  with 
the  method  of  alternating  directions.  This  approach  involves 
//steps  with  the/th  step  requiring  the  solution  of  theequation 


[*5  +PJL  +  K2  k2)]^- 

dr 

—  lktflj.„  (L  +  K  ~  —  k;,)U.  (14) 

We  solve  Eq.  (14)  by  first  discretizing  depth  dependence 
with  Galerkin’s  method  as  described  in  the  Appendix.  The 
resulting  system  is  then  solved  with  Crank-Nicolson  inte¬ 
gration. 

A  special  version  of  PE,  has  been  considered  for  appli¬ 
cation  in  shallow  water. l:  Since  sound-speed  variations  are 
very  small  in  shallow  water,  a  signal  trapped  in  a  shallow 
ocean  is  influenced  more  by  diffraction  than  refraction.  The 
shallow-water  version  of  PE,. 


II 

k„)U-  r,p\k 

2  ikttL  v 

(15) 

dr 

4k  r,  +  L 

obtained 

from  PE, 

by  assuming 

that 

|  (A" :  —  k  l )  U  |  4 \LU  |  in  the  water  column.  Equation  (15) 
has  the  same  refraction  term  as  Eq.  (9)  but  an  improved 
diffraction  term. 

To  illustrate  the  ability  of  PE„  to  handle  long-range  and 
very-wide-angle  propagation,  we  consider  a  waveguide  250 
m  thick  with  pressure-release  top  and  bottom  boundaries  in 
which  c=  1500  m/s.  A  25-Hz  point  source  is  placed  at 
z  =  25  m,  and  we  take  c„  =  1500  m/s.  The  eight  normal 
modes  for  this  problem  propagate  at  approximately  7, 1 4, 2 1 , 
29,  37,  46,  57,  and  74  uvg  from  horizontal.  The  PE„  solu¬ 
tions  (initialized  with  the  normal-mode  solution  at  r  =  0) 
are  compared  with  the  normal-mode  solution  in  Fig.  1.  We 
observe  that  the  PE„  solutions  break  down  very  rapidly  with 
r  for  small  n.  However,  the  PE,,  solution  is  nearly  perfect  at 
r  =  4  km. 

We  now  consider  a  realistic  example  that  illustrates  an 
application  of  PE„  for  low-frequency  underwater  acoustic 
propagation  in  deep  water.  In  the  water  column,  we  assume 
the  Munk  sound-speed  profile14 


ciz)  =  ctl, 


1  +/' 


+  exp 


(16) 


where  n  =  0.0071,  cL.h  =  1500  m/s,  zch  =  1 000  m.  and 
H  =  1200  m.  The  ocean  depth  is  5000  m.  In  the  sediment, 
c  =  1850  m/s,  p  =  1.5  g/cm',  and /?  =  0.5  dB//i.  A  10-Hz 
point  source  is  placed  at  z  =  200  m.  and  we  take  c„  =  1 500 
m/s.  The  homogeneous  half-space  field1'  is  used  to  initialize 
the  field  at  r  —  400  m.  The  Lloyd’s  mirror  beams  produced 
by  the  source  propagate  at  approximately  1 1, 34,  and  70  deg. 
PE,  should  accurately  account  for  the  1 1-  and  34-deg  beams 
for  well  beyond  r  =  20  km.  From  the  PE,  and  PES  solutions 
appearing  in  Fig.  2,  however,  we  observe  that  PE,  cannot 
handle  the  70-deg  beam,  which  is  partially  reflected  from  the 
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FIG.  1.  Transmission  loss  at  z  =  25  m  for  a  25-Hz  source  in  a  waveguide 
with  perfectly  reflecting  boundaries.  The  dashed  curves  are  the  normal 
mode  solution.  The  solid  curves  are  the  PE„  solutions  for  (a)  n  =  2,  (b) 
n  =  4,  and  (c)  n  =  6. 


FIG.  2  Transmission  loss  at  z  =  20()  in  for  a  10-Hz  source  in  deep  water 
The  dashed  curve  is  the  PE,  solution  The  solid  curve  is  the  PE,  solution 


II.  THE  HIGHER-ORDER  TDPE 

Just  as  refraction  is  less  important  than  diffraction  for  a 
signal  trapped  in  shallow  water,  the  effects  of  attenuation 
and  dispersion  are  less  important  than  refraction  and  dif¬ 
fraction  for  must  problems.  Thus  we  do  not  derive  correc¬ 
tions  for  the  attenuation/dispersion  operator  in  the  TDPE  of 
Ref.  1 2.  We  assume  that  K  is  real  ( no  attenuation )  and  that  c 
is  independent  of  w  ( no  dispersion )  in  the  analysis  and  write 
Eq.  (13)  as 

—  -  /c;  X  a,.z^+P-*L  a  (17) 

dr  jt'i  Vj,„(o  +  8jnL 


where 


>,.=^(4-4). 

c0  \c  cl) 


n.  =  4  +  6/-(4~4)' 

cl  \C  Co) 


T o  obtain  a  higher-order  TDPE  that  is  easy  to  solve  numeri¬ 
cally,  we  write  Eq.  (17)  in  the  form 


ocean  bottom  and  makes  a  significant  contribution  to  the 
field  for  5  km  <  r  <  1  5  km. 

Since  the  parabolic  approximation  is  based  on  an  expan¬ 
sion  about  a  reference  wavenumber,  it  is  not  obvious  that  this 
approximation  cun  be  generalized  to  elastic  media  in  which 
waves  of  different  speeds  exist.  To  investigate  the  possibility 
of  generalizing  PE„  to  handle  solid  ocean  bottoms,  we  con¬ 
sider  a  problem  for  which  c  —  1 500  m/s  in  the  water  and  the 
ocean  depth  is  200  m.  In  the  fluid  ocean  bottom,  c  =  1700 
m/s,  p  —  1.5  g/cm\  ..’H  0  -  0  5  dfl//  A  W.H7  point 
source  is  placed  at  z  =  25  m.  The  PE,  solution  for  c„  =  1 500 
m/s  and  the  PE  ,s  solution  fore,,  =  300  m/s  appear  in  Fig.  3. 
The  agreement  of  the  solutions  suggests  that  a  higher-order 
elastic  PE  based  on  the  Fade  series  would  handle  both 
compressional  and  shear  waves  simultaneously. 


FIG  V  Transmission  loss  al  z  25  m  for  a  50-Hz  source  in  shallow  w  ater  — — ' 

The  dashed  curve  is  the  PE.  solution  for  c„  1 500  m/s.  The  solid  curve  is  ' 
the  PI  1 .  solution  fore,,  200  m/s.  “ 
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1 


dU  .  " 

- =  no  > 

or 


v  ( °^-  ^  P±" _ (22) 

“i\ yiniL>:  +  5;„L  j 


We  define 


u(r,z,t) 


j  U(r,z,(t>)e 


xp(  —  uot)d(o 


and  invert  the  Fourier  transform  in  Eq.  (22)  to  obtain  the 
higher-order  TDPE: 

du  y  (a,„  |  (24) 

dr  ,  -  Y,.„  (d'/dt 1 )  +  8inL)  dt 

which  we  refer  to  as  TDPE„ . 

The  alternating  directions  solution  of  TDPE„  requires 
numerical  solutions  for  each  of  the  following  n  +  1  equa- 


Since  Eq.  (25)  is  similar  to  the  refraction  term  of  the  shal¬ 
low-water  TDPE,  and  Eq.  (26)  is  similar  to  the  diffraction 
term  of  the  shallow-water  TDPE,  the  numerical  solutions 
developed  in  Ref.  12  can  be  modified  slightly  to  obtain  the 
time-domain  solution  of  Eq.  (22).  In  contrast,  the  time-do¬ 
main  solution  of  Eq.  ( 17)  requires  the  numerical  solution  of 
n  equations  similar  to  Eq.  ( 26 )  as  well  as  n  third-order  equa¬ 
tions  that  are  much  more  complicated  than  Eq.  (25). 

As  in  Ref.  12,  the  source  function  f(t)  is  assumed  to 
have  compact  support,  and  a  time  window  /,</</,  that 
contains  the  signal  at  all  times  is  chosen.  The  boundary  con¬ 
dition  u  =  0  is  imposed  at  the  pressure-release  surface,  deep 
within  the  sediment  at  z  =  zM  from  which  no  energy  returns 
to  the  water  column  due  to  attenuation,  and  after  the  signal 
has  passed  the  observer  at  t  —  t2.  The  boundary  conditions 
u  =  du/dt  =  0  are  imposed  before  the  signal  is  detected  at 
1  =  1,.  Equation  (25)  is  a  first-order  hyperbolic  equation 
that  can  be  solved  with  the  Lax-Wendroff  scheme. Ih  Galer- 
kin's  method  is  used  to  discretize  depth  dependence  in  Eq. 

( 26 )  as  described  in  the  Appendix.  The  resulting  equation  is 
then  solved  with  Crank-Nicolson  integration  in  r  using  cen¬ 
tered  differences  in  t  while  sweeping  from  /  =  /,  to  /  =  r2- 

To  demonstrate  the  ability  of  TDPE,,  to  handle  very- 
wide-angle  propagation,  we  consider  a  waveguide  of  thick¬ 
ness  300  m  with  pressure-release  top  and  bottom  boundaries 
in  which  c—  1500  m/s.  The  Gaussian  source 
f(t)  -■  exp[  —  ( v/):  ]  is  placed  at  z  =  25  m,  where  v  =  150 
s  '.  The  image  solution,  which  is  exact,  is  used  to  initialize 
the  field  at  r  =  200  m,  and  we  take  c„  —  1500  m/s.  The 
TDPE,,  TDPE.,  and  TDPE,  solutions  are  compared  with 
the  imauc  solution  in  Fig.  Each  of  the  solutions  is  very 
accurate  for  the  first  ariivals,  which  propagate  at  small  an¬ 
gles  However,  the  agreement  imptoves  with  n  tor  the  later 
arrivals,  which  propagate  at  larger  angles. 

In  past  studies  of  the  TDPE.  a  stability  condition  for  the 
numerical  solution  of  the  refraction  operator  has  been  dis¬ 
cussed.  However,  the  numerical  solution  of  the  diffraction 
operator  appeared  to  be  unconditionally  stable  based  on  nu¬ 
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FIG.  4.  Time  scries  at  r  -  MX)  m  and  z  2(X)  m  for  a  Gaussian  pulse  in  a 
waveguide  with  perfectly  reflecting  boundaries  The  dashed  curves  are  the 
image  solution.  The  solid  curves  are  the  (a)  TDPE,.  (b)  TDPE-,  and  (c) 
TDPE,  solutions. 


merical  results  While  performing  the  calculations  for  the 
previous  example,  however,  we  observed  a  new  stability  con¬ 
dition  involving  the  grid  spacings  Az  and  A/.  In  a  homoge¬ 
neous  medium,  the  numerical  solution  of  the  diffraction  op¬ 
erator  is  unstable  for  Az  =  c„A/  and  n  >  1 .  The  solution 
appears  to  be  stable  for  all  n  if  Az  >  Ac,, A/,  where  numerical 
experiments  give  A  s  1.4. 

To  demonstrate  the  ability  of  TDPE,,  to  handle  large 
variations  in  sound  speed,  we  consider  an  ocean  of  depth  400 
m  in  which  c  increases  linearly  from  1500  m/s  at  z  -  0  to 
1600  m/s  at  z  =  400  m.  In  the  sediment,  c—  17(X)  m/s. 
p  —  1.5  g/cm\  and  (3  0.5  dB/x.  The  Gaussian  source 

function  with  v  =  150s  1  is  placed  at  z  =  50  m.  and  we  take 
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c0  =  1 500  m/s.  The  half-space  field  is  used  as  an  initial  con¬ 
dition  at  r=  100  m.  The  plane-wave  loss  operator  of  Ref.  12 
is  used  to  model  attenuation.  However,  we  have  found  that 
greater  accuracy  is  obtained  by  using  c  rather  than  c„  in  the 
loss  operator.  Sediment  dispersion  is  neglected.  The  re¬ 
sponse  to  /  is  convolved  as  in  Ref.  1 2  to  obtain  the  response 
to  a  50-Hz  time-harmonic  source.  Transmission  loss  for  the 
TDPE,,  PE,,  narrow-angle  PE,  and  shallow-water  PE,  solu¬ 
tions  appears  in  Fig.  5.  The  narrow-angle  PE  solution  has 


0  2  4  6  8  10 

(a)  Range  (km) 


I  I ( >  ?.  I  ransmissmn  loss  at  r  390  m  for  a  ?{)-Uz  source  in  a  refracting 
ocean  Hie  dashed  curves  are  the  I’F,,  solution.  The  solid  curves  aretlte  (a) 
narrow -angle  Tl  I  b)  shallow-water  I’t:,.  and  (e)  I  DPT,  solutions. 
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large  phase  errors.  The  shallow- water  PE,  solution  is  better, 
but  it  too  has  a  large  error.  The  excellent  agreement  of  the 
TDPE,  and  PE,  solution^  demonstrates  the  ability  of 
TDPE,  to  accurately  handle  pulse  propagation  in  deep  wa¬ 
ter  and  shows  that  the  plane-wave  loss  operator  is  accurate 
for  this  problem. 

III.  CONCLUSIONS 

The  phase  error  problem  of  the  PE  model  has  been  com¬ 
pletely  eliminated.  The  higher-order  PE  based  on  a  Pade 
series  handles  problems  involving  very-wide-angle  propaga¬ 
tion,  propagation  in  the  nearfield,  propagation  out  to  long 
range,  and  propagation  in  domains  in  which  sound-speed 
variations  are  large.  Numerical  results  suggest  that  the  ap¬ 
proach  might  be  useful  for  elastic  wave  propagation  for 
problems  involving  a  superposition  of  eompressional  and 
shear  waves.  The  higher-order  PE  was  solved  in  the  time 
domain,  and  a  new  stability  condition  was  observed. 

ACKNOWLEDGMENTS 

This  work  has  been  approved  for  public  release  and  was 
supported  by  the  Office  of  Naval  Research  and  the  Naval 
Ocean  Research  and  Development  Activity  (Contribution 
No.  221:021:89). 

APPENDIX:  DEPTH  DISCRETIZATION  WITH 
GALERKIN’S  METHOD 

Galerkin’s  method  is  effective  for  discretizing  the  rela¬ 
tively  complicated  depth  operators  in  Eqs.  ( 14)  and  (26), 
which  have  coefficient  functions  that  may  be  discontinuous 
and  involve  derivatives.  The  depth  grid  points  are  defined  as 
z,  —  i&z.  The  basis  functions  (z)  vanish  for  jz  —  z,  j  >  Az, 
increase  linearly  from  0  to  1  over  z, .  ,  <  z  <  z, ,  and  decrease 
from  1  toOoverz,  <z<z,  +  , .  Wedefine  u,  ( r,t )  =  u(r,z,,t) 
as  well  as  ©,  =  0(z, )  and  <t>,  =  4>(z, )  for  arbitrary  func¬ 
tions  ©  and  <J>.  The  basis  functions  provide  the  approxima¬ 
tions 


u(r,z,r)s^M,(/-,/)'P,(z), 

i 

(Al) 

< K:)sp,1t,(z), 

l 

( A2) 

0(z)  ==£©,'!',  (z). 

i 

(A3) 

The  following  orthogonality  condition  is 
all  /: 

required  to  hold  for 

("  ;./ >  h  ° 

( A4) 

Equation  (26)  is  discretized  by  substituting  Eq.  ( A 1 )  for  u 
and  Eqs.  (A2)  and  (A3)  for  (he  coefficient  functions  into 
Eq.  (A4).  The  following  approximations  are  obtained  for 
the  depth  operators: 
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©,  i  +  60,  +  0,  +  i 

'  "I  — - 


( A5) 


„  ,  ®.-.  +®r 

Ou  .  ,  s - u, 


d  :u 


12 

0, 


dr 
dQdu_ 
dz  dz 


20,  0, 

u,  H - ;■  u, 


,  0„,  +©, 

+ - — - «,  + 1 . 


(Az)-  (Az)J  '  (Az) 

(<t>,  ,  +  2<t>, )  (0,  ,  -0,) 

- - ; -  U,  , 

6(Az)- 

<t>,  .  (0,  -  0,  ,)  + 20,(20,  -0,  +0^,(0, -0,,,) 

+ - — — ; - u, 


6(A z)2 


(0,.  ,  +20,  )(©,,,  -0,) 

H - - Ut 

6  ( Az ) 2 


( A6) 


(A7) 


These  approximations  are  also  used  to  discretize  depth  in 
Eq.  (14).  As  in  Ref.  12,  0  corresponds  to  leg(  p)  in  Eq. 
(A7). 
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